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Abstract 

Problems  on  the  identification  of  two-dimensional  spatial  domains  arising  in  the 
detection  and  characterization  of  material  damage  are  considered.  For  electromag¬ 
netic  nondestructive  evaluation  systems,  observations  of  the  magnetic  flux  from  the 
front  surface  are  used  in  a  output  least-square  approach.  Parameter  estimation  tech¬ 
niques  based  on  the  method  of  mappings  are  discussed  and  approximation  schemes 
are  developed  applying  a  finite-element  Galerkin  approach.  Theoretical  convergence 
results  for  computational  techniques  are  given  and  results  are  applied  to  numerical 
experiments  to  demonstrate  the  efficacy  of  the  proposed  schemes. 


1  Introduction 

Detection  and  characterization  of  small  cracks  and  corrosion  embedded  in  the  structures 
of  aircraft  are  critical  issues  in  the  maintenance  of  aging  aircraft.  For  instance,  carbon 
fiber  reinforced  plastics  (CFRP)  have  been  widely  used  in  the  advanced  aircraft  and,  as 
a  result,  demand  has  grown  for  assessing  the  structural  integrity  of  structures  made  from 
those  materials.  An  important  effort  on  such  problems  entails  quantitative  nondestructive 
evaluation  methods  in  SQUID-based  NDE  system  [1],  It  is  well  known  that  SQUIDs  (su¬ 
perconducting  quantum  interference  devices  )  have  very  high  sensitivity  for  magnetic  flux 
measurements.  Due  to  the  highly  sensitive  magnetic  flux  measurements,  inverse  analyses 
on  electromagnetic  problems  are  effectively  used  for  detecting  and  characterizing  cracks 
and  delamination.  Initial  efforts  on  such  inverse  problems  include  nondestructive  testing 
under  SQUID  measurements  using  the  direct  current  method  [2],  However,  the  detections 
of  non-through  crack  and  corrosion- like  damage  by  direct  current  flows  are  insufficient  be¬ 
cause  of  the  lack  of  information  for  deep-lying  effects  from  vertical  component  of  magnetic 
flux  measurements.  Since  a  skin  depth  of  metal  varies  according  to  frequency  of  flowing 
current,  it  is  possible  to  detect  the  deep-lying  flows  by  controlling  the  frequency.  To 
this  end,  SQUID  based  nondestructive  evaluation  (NDE)  systems  using  injection  current 
methods  have  been  recently  developed  [3,  4],  In  this  paper,  we  propose  a  computational 
method  for  recovering  corrosion- like  damage  with  magnetic  flux  density  data  from  the 
SQUID  based  NDE  system  with  alternating  current  force.  The  idea  of  the  method  of 
mapping  techniques  ([5,  6])  are  effectively  used  in  our  inversion  procedures.  Figure  1 
illustrates  the  overall  configuration  of  SQUID  based  NDE  system  proposed  here.  In  this 
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inspection  procedure,  alternating  current  forces  with  multiple  angular  frequencies  are  ap¬ 
plied  to  the  sample  material  inspected.  A  damage  corrupts  the  current  flows  inside  the 
conductor  and  we  can  detect  this  disturbance  from  the  magnetic  flux  measurements  by 
the  SQUID.  The  task  we  consider  here  is  to  identify  the  geometrical  shape  of  the  damage 
from  input  and  output  data.  Throughout  the  paper,  we  assume  that  the  damage  is  lo¬ 
cated  sufficiently  far  from  both  sides  of  the  sample  and  that  the  direction  of  current  flow  is 
uniform  with  respect  to  the  length  direction.  In  Section  2,  we  consider  a  two-dimensional 
mathematical  model  for  inspections  based  on  a  SQUID-based  NDE  system.  In  Section  3, 
a  direct  problem  is  formulated  in  variational  form  in  an  appropriate  Hilbert  space.  In 
Section  4,  the  inverse  problem  is  discussed  in  the  context  of  a  nonlinear  output  least 
square  problem.  The  class  of  admissible  parameters  are  given  and  the  existence  of  solu¬ 
tions  is  demonstrated  using  the  idea  of  the  method  of  mappings.  Section  5  is  devoted  to 
theoretical  convergence  for  the  proposed  computational  method.  Numerical  experiments 
are  summarized  in  Section  6  to  demonstrate  the  efficacy  of  the  proposed  method.  Related 
efforts  on  similar  problems  are  given  in  [7,  8,  9]  where  eddy  current  based  techniques  are 
used  and  the  focus  is  model  reduction  techniques  are  opposed  to  the  use  of  the  method 
of  maps  employed  here. 

2  Mathematical  Model  of  Inspection  Procedures 

Our  forward  analysis  is  considered  on  a  cross  section  of  the  conducting  material  sample, 
he., 

Gc  =  {  x  =  (x1,x2)  |  \xx |  <  dci  0  <  x2  <  h0  } 
as  illustrated  in  Fig.  2.  Assuming  that  the  material  is  of  infinite  extent  in  the  x3  direction, 
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Figure  2:  Inspection  Area 


our  problem  is  defined  on  an  appropriate  computational  “window”  given  by 


G  =  {  x  |  |xi|  <  d,  |a:2 1  <  h  (ho  <  h  <  oo)  } . 


A  set  of  alternating  current  sources  J .-L  with  frequencies  u is  applied  to  the  sample  material 


J i(t,x)  =  (0,0,  Js(x))T  cos(cjit). 


The  set  of  applied  frequencies  u 
depth  condition, 


can  be  determined  in  accordance  with  the  skin 


(1) 


where  /i  denotes  the  magnetic  permeability  and  where  a  denotes  the  electrical  conductivity 
of  the  sample.  Given  that  skin  depth  of  a  metal  changes  according  to  frequency  of 
flowing  current,  it  is  possible  to  detect  a  configuration  of  deep-lying  defects  by  varying 
the  frequency.  For  the  output,  measurements  can  detect  the  gain  margin  and  phase 
change  of  the  voltage  using  a  LCZ  meter.  This  means  that  the  SQUID  measurements  can 
be  represented  as 

B (£,  x)  =  B(x)  cos(c jt  +  0(x)). 


Thus  it  is  natural  and  convenient  that  the  state  variables  be  described  by  complex  phasor 
representations.  The  phasor  form  of  Maxwell’s  equations  is  given  by 


V-B  =  0, 


(2) 
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V-E  = 

0, 

(3) 

V  x  E  = 

~ju  B, 

(4) 

V  x  H  = 

J. 

(5) 

Following  a  standard  approach,  we  introduce  the  magnetic  vector  potential  A  = 

(AXl,  AX2,  Ax 3)  defined  by  B  =  V  x  A  and  the  cross  product  of  Eq.  (  4  )  can  be  replaced 
by  representing  E  +  juA  as  the  gradient  of  a  scalar  electrical  potential  0,  i.  e., 

E  =  —ju  A  —  V0. 


In  conjunction  with  Ohm’s  law  J  =  a E  and  the  constitutive  law  B  =  pH,  the  system 
state  A  is  governed  by 


1 

V  x  -V  x  A 

h 

V  •  ( juA  +  V0) 


—j  luctA  —  aVc/) 

0. 


(6) 

(7) 


We  assume  here  that  the  sample  is  of  nonmagnetic  material,  so  that  the  magnetic  per¬ 
meability  p  is  equal  to  that  of  air,  i.e.,  p  =  p0.  In  our  formulation,  the  magnetic  vector 
potential  becomes  A  =  (0,  0,  A3).  Therefore,  the  equation  for  the  component  of  magnetic 
vector  potential  A3  =  A  is  simply  rewritten  as 

1  _2  ,1  A  <90 

- V2A  +  jXcVuA  =  -Xccr^— 

Po  0x2, 


where  Xc  denotes  the  characteristic  function  of  the  conducting  region  Gc.  Taking  into 
account  that  the  right  side  of  the  above  equation  is  composed  of  the  source  current 
density,  it  follows  that 


=  —a 


<90 

<9X3 


in  Gr 


Consequently  we  can  formulate  the  two  dimensional  forward  problem  in  terms  of  the 
equation: 


1 

ho 


f  <92A(xi,x2)  <92A(xi,x2)\ 

y  <9x2  <9x2  ) 


+  jXc(?uA(x  i,x2) 


XcJs{x  1,X2). 


(8) 


Since  the  boundaries  in  the  x2  directions  are  assumed  to  be  sufficiently  far  from  the 
sample,  we  assign  a  Dirichlet  boundary  condition  on  this  part  of  the  boundary,  i.  e., 


A  =  0  on  Se> 


(9) 


where 


Sd  —  {  x  |  |xi|  <  d,  \x2\  =  h  }. 

We  also  consider  the  boundaries  in  the  X\  directions  to  be  sufficiently  far  from  the  damage 
so  that  Neumann  boundary  conditions  in  the  X\  direction  can  be  set  as: 


<9A 
<9  n 


0  on  Sn 
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(10) 
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Figure  3:  Setting  of  the  boundary  condition  on  the  window 


where 


SN  =  {  x  |  |xi|  =  d,  \x2\  <  h  }. 


Figure  3  denotes  the  boundary  conditions  on  the  computational  window.  In  contrast 
to  other  magnetic  sensors,  SQUID  measurements  can  directly  detect  the  magnetic  flux 
densities  that  are  independent  of  the  applied  frequencies.  Suppose  that  the  scanning 
strategy  is  given  on  the  line 


/  =  {  x  G  G  |  \xi\  <  r  (dc  <  r),  x2  =  —hf  }, 
as  shown  in  Fig.  4.  Thus  the  observation  mechanism  involves  data 

BA 

B2(x.p)  =  —  — — (xp)  for  xp  G  I,  p  =  1,  2,  •  •  • ,  m. 

uX\ 

and  the  observations  are  given  by 

z  =  HA  eCm 


where  7i  denotes  the  bounded  operator  given  by 


mr  = 


d(- 


\Op\  Jop  dxi 
and  Op  denotes  the  sub-interval  of  /  defined  by 


dx-i 


X2=—hf 


(11) 


(12) 


(13) 


Op  C  (p  —  1,  2,  •  •  ■ ,  m). 

and  \Op\  denotes  the  length  of  the  sub- interval  Op. 


(14) 
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Figure  4:  Scanning  procedures 

3  Weak  Formulation 

Let  q  be  a  constant  vector  which  characterizes  an  unknown  damage  shape  where  we 
assume 

[H-0]  The  admissible  set  Q  of  damage  parameters  is  a  compact  subset  of  RM  . 

We  introduce  the  “suspect  region"  Gs  C  G  which  covers  the  suspected  but  unknown 
damage.  More  precisely,  as  depicted  in  Fig  5,  the  domain  Gs  is  described  by 

Gs  =  {  x  =  (xi,  x2)  |  |xi|  <  c,  0  <  x2  <  h  } . 

For  convenience  of  discussions,  the  suspect  region  Gs  is  decomposed  into  two  subre¬ 
gions  which  correspond  to  the  air  region  G*  and  the  conducting  region  G2S  and  these 
are  parametrized  by  q.  Thus  the  “window”  G  is  decomposed  into  three  sub-regions  Go, 
G]  and  G2S  as  follows: 


Gs  —  G*(q)  U  Gy(q),  G0  =  G  -  Gs. 

We  assume  the  boundary  of  G  is  decomposed  into 
S°N  =  {x  |  |xi|  =  d,  \x2\  <h} 

Sp  =  {x  |  c  <  |xi|  <  d,  x2  =  h  }  U  {x  |  |xi|  <  d,  x2  =  —h  } 
Sp  =  {x  |  |xi|  <  c,  x2  =  h  } 
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X2 


and  that  the  interface  among  sub-regions  are  defined  by 
Sj1  =  {x  |  |xi|  =  c,  h0  <  x2  <  h  } 

S?  =  {x  [  \x\\=  C,  0  <  X2  <  ho  }  U  {x  |  |xi|  <  C,  X2  =  0  } 

and  S}2  is  defined  by  the  interface  between  G](q)  and  Gpq).  The  form  of  S}2  will  be 
described  in  the  later  discussions.  The  system  state  is  then  denoted  by 

A  :=  A0  in  G0  A  :  =  A{  in  G*  (i  =  1,  2) 

where 

A-Af+jAl  (z  =  0,1,2). 

In  the  sequel,  we  suppose  that  the  electric  conductivity  of  the  conducting  material  sample 
has  constant  value  ac,  i.e. 

®  XcVc- 

Thus  the  subsystem  on  Go  is  described  by 

— -V2A*  -  Xc<?cMjA{}  =  XcJs  (15) 

Ho 

— -V2Aq  +  XcVcUjAq  =  0  in  Go  (16) 

Ho 
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with  the  boundary  condition 


K  =  K  =  « 


BAS-  BM 


=  0 


on  Sp 
on  S°N. 


dn  dn 

The  subsystems  defined  on  G*(q)  and  G2(q)  are  described  respectively  by 

V2Af  =  0 

W'2A[  =  0  in  Gj(q) 

with  the  boundary  condition 

A?  =  A[  =  0  on  Sjj, 

and 

1 


ho 

1 

ho 


-V2Llf  —  o,jjjA =  Js 


V2A2  +  ocijjA2  =  0  in  G2(q). 


We  also  need  the  interface  conditions  between  G0,  G^(q),  and  G2(q), 
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5J 


U2/ 


G0UGl(q)UGl(q) 


|V<h|  dx\dx2. 


P  :=  (<hRV)  ehgx  Vq. 

Then,  for  any  tp,  /  €  Vq  x  Vq,  we  define  the  bilinear  form  on  Vq  x  Vq 

A  ft 

ho  J  JGq\JG\  (q)UG|  (q) 

XcVcVIi,Rdxidx2 


«( q)(&^)  := 


(\7pR  •  Vi/  +  Vp1  •  V/1)  dXldx2 


—U) 


GoUGi(q) 


+uj  1 1  Xc®cP  4>  dx\dx2 

/G0UG2(q) 


(17) 

(18) 

(19) 

(20) 

(21) 

(22) 

(23) 

(24) 

(25) 

(26) 

(27) 

(28) 
(29) 


(30) 


and  the  linear  form  on  V„ 


L( q)(<p)  :=  //  XcJs^Rdx1dx2. 

J  JGqUG2  (q) 


(31) 


Lemma  1:  Suppose  that 

[H-l]:  The  applied  angular  frequency  u>  is  properly  chosen  such  that 

OJ  <  CUo  <  oo 

and  that 

XcJs  e  L2(G). 

Then,  for  every  q  G  Q,  there  exists  a  unique  solution  A(q)  =  ( AR ,  A1)  EVqxVq  which  is 
the  solution  of 

a(q)(A  p)  =  L(q)(p)  for  p  G  Vq  x  Vq.  (32) 

Moreover,  we  have 

ll^(cl)l|vrgxV9  <  Ki\\XcJs\\l2(G)  (33) 

where  K\  is  a  constant  independent  of  q. 

Proof:  For  arbitrary  p  e  Vq  x  Vq,  a( q)  is  Vq  x  V^-elliptic  with  constant  5  =  l//io  since 


a(  q)(<p,<p)  = 


WpR  '  +  Wp1  '  dxtdx2  =  d\\(pfVaXVa  (34) 


hO  ^^G0UGl(q)UG2(q) 

from  which  follows  the  coercivity.  For  the  boundedness,  it  can  be  easily  checked  that 


a{ q)(&^)  <7 


(35) 


The  linear  functional  L(q)  is  bounded  on  Vq  x  Vq,  i.e., 


lL(q)(^)l  = 

< 

< 


//  XcJsPRdxldx2  +  //  JspRdx1dx2 
JJgo  jJg2s{  q) 

/ /  XcJsPR  dxldx2  +  //  dxidx2 

JJgo  JJg2s(  q) 


Go 


|Xc</s|  dx\dx2  + 


G?(q) 


I JJ2  dx\dx2 


x||^||L2(GoUGi(q)uG2(q)) 


<  M 


for  p  £  Vq  X  Vq. 


Applying  the  Lax-Milgram  lemma,  for  each  q  e  Q,  we  hnd  there  exists  a  unique  solution 
A(q)  G  y,  x  Vq  satisfying  (32).  Similarly,  we  have 


^P(q)llv0xv0  <  a(q)(A^)  <  MP(q)llvgxvg 
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X2 


Figure  6:  The  unknown  damage  function 

from  which  follows 

M 

IKq)lkxV,  <  K±  where  K\  =  — . 

This  completes  the  proof. 

Remark:  In  [H-l],  the  upper  bound  of  angular  frequency  ujq  can  be  reasonably  set  for 
the  inspection.  More  precisely,  the  angular  frequency  uj  of  the  applied  current  force  must 
be  properly  chosen  by  the  skin  depth  condition  (1). 

4  Admissible  Class  and  Method  of  Mapping 

In  this  section,  we  restrict  the  corrosion-like  damage  so  that  it  can  be  represented  by  a 
simple  parametrized  function.  The  unknown  defect  shape  is  characterized  by  a  function 

x2—r(x  i ;  q)  for  |xi|  <  c  (36) 

with 

0  <  r(x i;  q)  <  h0  (  <  h  <  oo  )  and  r(— c;  q)  =  r(c;  q)  =  h0 

as  depicted  in  Fig.  6.  We  describe  by  Q  C  RM  an  admissible  set  of  possible  parameter 
values.  The  geometrical  structure  of  the  damage  is  restricted  by  the  following  hypotheses: 

[H-2]  For  each  q  G  Q,  r(-,  q)  G  W4(— c,  c)  . 
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Figure  7:  Refence  domains 

[H-3]  There  exists  a  constant  (3  such  that,  for  each  q  e  Q, 

0  <  (3  <  r(-,  q). 


[H-4]  There  exists  a  function  dist  :  Q  x  Q  — >  R1  with  dist( q,  q)  — >  0  as  |q  —  q|  — >•  0  such 
that 

lk(-,q)  -  ^(-,q) II  1,00  <  dist(q,q)  for  q,q€Q. 

The  subregions  Gi(q)  and  G2(q)  given  in  the  previous  section  can  now  be  explicitly 
rewritten  as 

G^(q)  =  {  x  |  |xi|  <  c,  r(x i,q)  <  x2  <  h  } 

and 

G^(q)  =  {  x  |  \x\  |  <  c,  0  <  x2  <  r(xi,q)  }  , 

respectively.  Associated  with  these  sub-domains,  we  define  the  reference  sub-domains  as 
shown  in  Fig.  7: 

=  {  x  |  l^i  |  <  c,  ho  <  x-i  <  h  } 

G2s  =  {  x  j  |xi|  <  c,  0  <  X2  <  ho  } . 

We  introduce  the  affine  mapping: 


x  =  T(q)x 


Ti(q)x  for  x  G  G] 
T2(q)x  for  x  G  G2S 


(37) 
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Figure  8:  Transformation  mapping  from  the  reference  domain  into  the  domain  with  dam¬ 
age 


where 


Ti(q) 
U  q) 


Xi  =  X\ 

x2  =  (h -r(xi,q))(x2  -  h)/(h  -  h0)  +  h 

X\  =  X\ 

x2  =  r(x  i;q)x2/h0 


These  mappings  result  in  the  following  identities  : 

Gi(q)=Ti(q)oGj  (*  =  1,2). 


(38) 


Figure  8  illustrates  the  transformation  from  the  reference  domain  into  the  domain 

with  damage  G*(q)  U  G^(q).  Let  (p  be  in  the  set 

V  =  {  (p  e  H\G0  U  G]  U  G])  |  <p  =  0  on  S°D  U  S]D  } 

endowed  with  the  norm 

ll^ll2  :=  If  _  .  \\/<p\2  dx1dx2. 

JJg0uGIug2s 

Let  us  set 


$  :=  G  V  x  V. 

Then,  for  any  <p,i/j  £V  x  V,  the  bilinear  form  a(q)  can  be  represented  on  V  X  V  as 
a(q)(|,|)  :=  {(VvR)-(^(q))"'(^Ti(cl))"1(^R) 

H'  0  {= i  J  V 
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(W)  •  (VTi(q))  f(vr4(q))  1  (V,0'r)|  det  |VTj(q)|  dx\dx2 

jJSgq{(?*R)  •  (W*)  +  (V^)  •  (y^)}dxl(ix2 


—00 


acipIripR  det  VT2(q)  dx\dx2 


G? 


+oj  1 1  _  acp  tp1  det  VT2(q)  dx\dx2 


-oo  il  Xc&c+  tpKdx\dx2  +  oo  //  Xc&c+  tp  dx\dx2. 

I  Go  JJGo 


The  linear  form  on  V  is  also  represented  by 

L(q)(<p)  :=  IJ  JsifR  det  VT2(q)  dx1dx2  +  JJ  XcJs+Rdxxdx2 


Noting  that 


V7K  q)  = 


/  9xi 

'  9xi  9x2 


9X2  9x2 

9xi  9x2 


we  obtain 


VTi(q)  = 
Similarly,  we  have 


/  1 
v  - 


(  1 


r'(xi;cj)(x2-h)  h—r(x  i;q) 

h—ho  h—ho 


VT2(q)  = 


r'(xi;q)x2  x(xi;q) 
ho 


ho 


det 


VTi(q) 


h  —  r(x i;  q) 


det 


VTx(q) 


h  —  ho 

Hence  the  explicit  form  of  Eq.  (39)  can  be  represented  as 


r(x i;  q) 

hn 


Kq)(M)  := 


i=l 


i  /  <9^  ch/)R  dip1  dip1 ' 
Ql  ^  1  <9xi  <9xi  cbq  <9x] 


.  /  <%5R  <%5R  dtpR  dip1  dip1  dip1  dtp1 ' 

+a2(q)  |  — — - b  — — xx - b  - b 


<95q  (9x2  (9x2  <9x]  dx i  <95i2  <95i2  <95q 

dxidx2 


+  *  (  )  (^ 


(39) 


(40) 


<95i2  (9x2  dx2  dx2 

dipR  dtpR  dip1  dtp1  dipR  dtpR  dip1  dtp1! 

po  JJg0  \  dx i  <9xi  <9xi  <9xi  <9ai2  cbc2  <95i2  <9ai2  J 

+  jj  ao(q)  {^-ipItpR  +  ipR tp 7)  dx\dx2 

+C9  Xc^c  (-PP1^11  +  ^tp1)  dxxdx2  (41) 
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where 


al(q) 

a2(q) 

«s(q) 

«i(q) 

a!(q) 

«s(q) 


r'(x i;  q)2(f2  -  h)2  +  (h  -  r(xi;  q))2 
Ho(h  -  h0)(h  -  r(x x;  q)) 
^(^i;q)(^2  -  h) 

/x0(/i  -  r(fi;q)) 
h  —  ho 

Ho(h  -  r(xi,q)) 
r'(x i;  q)2x|  +  r(xi;  q)2 
/x0h0r(xi;q) 
r;(gi;q)g2 
/i0r  (xi;q) 
ho 

AW(zi;q)’ 


and  where 


a0(q)  = 


a>crer(^i;q) 

ho 


From  (40),  the  transformed  linear  form  can  be  explicitly  rewritten  by 


£(q)(<l) 


Jsr(x i;  q) 
ho 


LpRdxidx2  + 


XcJs(fRdxidx2- 


(42) 


Lemma  2:  With  the  hypotheses  [H-0]  to  [H-4],  there  exist  positive  constants  cqA,/l2 
and  K:i  such  that,  for  qi,q2  €  Q,  the  bilinear  form  a(q)(-,-)  satishes  the  following  in¬ 
equalities  for  all  G  V  xV: 


«(q )(<P,<P)  >  ol 


VxV 


a(q)(frfy 

h(q)(<f,^)  -  a( q)(<p,$) 


<  k2 


VxV 


VxV 


<  K3dist( q,  q) 


VxV 


VxV 


(43) 

(44) 

(45) 


where 


dist(q,  q)  ^  0  as  |q  —  q|  — ^  0. 

Moreover,  a,  K2,  and  K3  can  be  chosen  as  constants  which  are  independent  of  the  pa¬ 
rameter  vector  q. 


Proof:  It  can  be  easily  shown  that,  with  Ci,C3  >  0,  any  quadratic  form  satishes 
Ci(^i)^  —  2c2^i^2  +  03(^2)  > — 3  2^  ^  min  { (ci)  \(c3)  (|£i|2  +  |^2 12)  - 
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From  (41),  the  associated  quadratic  form  for  each  suspect  region  Gls(i  =  1,2)  becomes 

{“l(q)Kf)2  +  2a- (q)KKf )  +  a-  (q)K2R)2} 

+  {“i(q)({i)2  +  2“2<q)te'?2)  +  «3(q)(?2)2} 

>  n'i(q)tt3<q)  -  (tt2<q))2  mln  {(a-(q))-1,  (a- (q))-1} 

x(|ef|2  +  |^|2  +  |^|2  +  |^|2)-  (46) 

We  easily  find  that 

ai(qK(q)  -  (<4(q))2  =  1  for  i  =  2- 

Under  the  hypotheses  [H-2]  and  [H-3],  we  admit  that 

Ho(h  -~  h0)(h  -  r(x pq))  ^  ho 


(«i(q))  1  = 


> 


r'(x i;  q)2(x2  -  h )2  +  ( h  -  r(x p  q))2  Mx  +  M2 


where 


-  Mo(h-r(x i;q)) 

(a3(q,) - jrnr„ —  -  ',0 


M\  =  snp  sup  |  r/(d:i ;  q)  | 3  M2  = 

qeQ  x\  e[— c,c] 


'  h-0 

h  —  tin 


>  1. 


This  implies  that,  for  G\ , 


min  {(ai(q))  \  (<4(q))  *} 

Similarly,  noting  that 

(a?(q))_1  = 


> 


ho 

Mi  +  M2 ' 


li0h0r(x  i;q) 


> 


ho^ 


(«a(q))  = 


r'(xi;  q)2X2  +  r(xi;  q)2  ho(Mi  +  1) 
ho'K5  i;q)  yU0/3 


we  obtain 


mm 


ho  ho 

{(a?(q))'1,(«3(q))'1} 


> 


/ip^  fnr  n 2 

ho(Mi  +  1) 


Choosing  the  constant  as 


a  =  mm 


ho 


Mo  P 


(47) 


(48) 


M\  +  M2  h0{M1  +  1)  /io  J 

which  is  independent  of  q,  we  obtain  the  coercivity  (43)  of  the  transformed  sesquilinear 
form.  To  prove  the  boundedness,  we  note  that,  under  the  hypotheses  [H-l]  to  [H-4],  it 
follows  that 


<?Sq  Qi 


a*  (q)  <  sup  sup  a) (q)  (=  M2)  for  i  —  1,  2\j  —  1,  2 


h  —  hn 


hn 


4(q)|  <  ^ 

|«o(q)|  <  cu0(jc  (=  m4). 


(=  M3)  for  i  —  1,2 
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Thus,  by  setting  as 


K2  =  max  (m2,  M3,  M4j  /j0  , 

we  obtain  the  boundedness  (44).  To  establish  the  continuity  property,  we  note  that,  for 
any  q,  q  G  Q, 


a(q)(&V0  -  a(q)(<p,$) 


dpR  dijjR  dip1  d'l/j1 


dx\  dx\ 

R  P)nI,R 


dx\  dx\ 


dpR  dtyR  |  dpR  d^R  |  dip1  d'd1  |  dp1  di^1 ' 


<  E  [L  {«i  (q)  -  °i  (q) } 

i=\  JJGl  L 

+  {a*2(q)  -«2(q)} 

+  (4(q)  -  4 (q)}  j  j  ditdi2 

jj^2  (ao(q)  -  o0(q)}  (-<E0R  +  dzicS 


<95;  i  dx2 
d(pR  d^11 


dx2  dx\  dx\  dx2  dx2  dx4 
dip1  dtjj1 ' 


With  the  hypotheses  [H-l]  to  [H-4],  it  can  be  argned  that  a*(q)(i  =  1,2]  j  =  1,2,3)  are 
continuous  in  L°°(G*).  We  can  thus  infer  the  continuity  of  the  bilinear  form  (41)  with 
respect  to  q  G  Q.  The  proof  has  been  completed. 


Lemma  3:  Let 

!(q)  =  Ao  T(q) 

be  the  transformed  system  state.  Then,  with  the  hypotheses  [H-0]  to  [H-4],  for  every 
q  G  Q,  there  exists  a  unique  solution 

l(q)  ebxh  (49) 


in  the  sense  that 

a(q)(i(q),<!)  =  L( q)(p)  for  p  e  V  x  V.  (50) 

Moreover  the  solution  kl(q)  in  the  system  on  G  =  Go  U  G7  U  G2S  is  bounded  in  V  x  V 
uniformly  in  q  e  (). 


Proof:  From  Lemma-2  and  from  (30),  (31),  (41),  and  (42),  we  have  that,  for  ip  G  Vq  x  Vq 
and  <p  G  V  x  V , 


a 


2 

VxV 


<  |a(q)(<^,  p)\  <  K2 


2 

VxV 


and 

=  a{q)(p,p)- 


Taking  into  account  that  the  constants  a  and  K4  are  independent  of  the  unknown  param¬ 
eter  q,  this  implies  that  the  Vq  x  V^-norm  is  equivalent  to  the  norm  of  V  x  V  uniformly 
in  q  e  Q-  Since,  under  the  hypothesis  (H-l), 


£(q)(<l) 


<  K 4 


L2(G)xL2{G)  ’ 
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there  exists  the  solution  A  in  V  x  V  in  the  sense  of  (50).  The  proof  has  been  completed. 

Lemma  4:  With  the  hypotheses  [H-0]  to  [H-4],  q  — >  A(q)  is  continuous  from  Q  to 
V  xV. 

Proof:  Let  qfc  — »  q  in  Q  and  let  A(qfc),  A(q)  be  the  corresponding  solutions  of  (50). 
Then,  from  (50),  we  have 

a(qfc)(i(qfe),<^)  -  a(q)(!(q),  <p)  =  [L( qk)  -  L( q)](£)  for  \/(p  G  V  X  V. 

Then  we  may  rewrite 

a(qfc)(I(qfc)  -  A(q),<p)  +  [a(qfc)  -  a(q)](l(q),  <p)  =  [L( qk)  -  L{ q)](^)  for  \/ip  E  V  x  V. 

Choosing  as  (p  =  A(qfc)  —  A(q)  in  V  x  V  and  using  (43)  and  (45)  in  Lemma-2,  we  find 
that 


«l^(q  )  -^(q)lvxv  -  dist( q  ,q)  (IxjsIl2(g)  +  l^(q)lvxvf  l^(q')  -^(q)lvxv 

This  yields  the  desired  continuity,  that  is, 


a|A(q')  -  v4(q)|pxy  <  dist( q  ,  q)  +  l^(q)lyxy 

The  proof  has  been  completed. 

The  observation  for  the  transformed  state  A  can  be  rewritten  by 

z(q-u)  =  HA(q;u) 

where  H  :  V  x  V  — >  Rm  x  Rm  is  given  by 

Iff  d(pR  dip1' 


0  as  q^ 


Hp> 


\Op\  Jop  l  dxi  ’  dxi 


dx  i 


(P=  1,2 


X2 =—hf 


q 


(51) 


Figure  9  depicts  the  procedures  for  gathering  the  measurement  data  from  SQUIDs. 

Lemma  5:  With  the  hypotheses  [H-0]  to  [H-4],  the  mapping  q  — >  z(q)  is  continuous 
from  Q  to  Rm  x  Rm. 

Proof:  Taking  into  account  that  the  operator  H  belongs  to  C{V  x  V ;  Rm  x  Rm )  for  each 
q  e  Q,  we  see  that  the  above  statement  follows  directly  from  Lemma-4. 


The  parameter  estimation  problem  is  formulated  as  follows: 


17 


X2 

A 


-  d 


h 

~1 

Gs 

ho 

-  r 

C'2 

Gs 

r 

-  dc  - c  O 

1  !  1  1 

±  c 

1  1  1  1 

\/^of 

Zi 

1  ^ 

'  nf  Scanning  line  I 

-h 


xi 


Figure  9:  Data  acquisition  by  SQUIDs 


Shape  Identification  Problem  SIDP: 


Given  the  observed  data  {^(u^)}^  from  applying  the  set  of  alternating  injection  currents 
{Js  cos(cj,d)}^lrj,  then  hnd  the  optimal  q  =  q*  which  minimizes  the  fit-to-data  functional 


1  Mo 

F(q)  =  9  51 rt  r)  ~  zd(^i) 

2  i= i 


(52) 


with  respect  to  q  G  Q  where  Q  is  a  compact  set  of  RM . 


Theorem  1:  With  the  hypotheses  [H-0]  -  [H-4],  the  problem  SIDP  has  at  least  one 
solution  q*  G  Q. 


This  follows  immediately  from  the  continuity  properties  given  in  Lemma  5  and  the 
compactness  of  Q. 


5  Computational  Method  and  Convergence  Analysis 


The  computational  scheme  is  based  on  the  use  of  a  finite  Galerkin  approach  to  construct 
a  sequence  of  finite  dimensional  approximating  identification  problems.  To  approximate 
SIDP.  we  choose  a  sequence  of  finite  dimensional  subspace  HN  x  HN  c  V  x  V  such  that 


nV-<^ 


V  xV 


0  as  N  — *  oo  for  (p  e  V  x  V 


(53) 


18 


where  IIA  is  the  orthogonal  projection  of 

L2(G0  UGjU  G2)  x  I/2 (G0  UG^IJ  G2)  onto  HN  x  HN  (54) 

-*N 

The  approximating  system  is  dehned  for  A  G  HN  x  HN  by 

a(q)(A  ,(f  )  =  L(q)((p  )  for  ip  G  HN  x  HN .  (55) 

The  observation  output  for  the  approximating  system  can  be  represented  as 

_  -zN 

zN (q; u)  —  HA  (q;o;).  (56) 

Thus  the  computational  method  is  formulated  as  follows: 

Approximate  Shape  Identification  Problem  (ASIDP)iV: 

Find  G  Q  which  minimizes 

1  m„  2 

^(q)  =  X  Z)  |^(q;  ui)  ~  zd(^i)  |  (57) 

i— 1 

Lemma  6:  Let  qA  ^  q  E  Q.  Then 

-*N 

A  (qN)  - >  i(q)  G  V  X  V.  (58) 


Proof:  We  have 

->N  ->  -> N  -  - 

A  (qN)~A(d)  <  A  (qN)-nNA(q)  +  UNA{q)  -  A(q) 

V  xV  VxV 

From  (53),  it  suffices  to  prove 

- N 

A  (qN)  -  UNA(q)  ^  0  for  q^  ->  q  G  Q  as  N  ->  oo. 

VxV 

We  obtain 

-.iv  ^ 

~a(qN)(A  (qN)^N)-d(q)(A(q)^N)^[L(qN)-L(q)}^N)  for  t/,N  G  HN  x  HN . 

Furthermore  we  have 

->N  -  - 

d(qN)(A  (qN)  -  UNA(q),^N)  +  a(q")(n"i(q)  -  A(q)^N) 

+  [5^)  -  a(q)](A(q)^N)  =  [L{ qN)  -  L(q)](^N). 
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_  -z 

Choosing  i/)N  =  AN  =  A  (qN)  —  II "^(q),  we  find  that 
«||AiL*  ^  K*  n^q)-^)  .  JAN\\yxy 

VxV  VxV 

+dist( qN,q)  A(q)  _  _  \\AN \\yxy  +  Mdist(qN ,  q)\\AN \\yxy. 
vxv 

Consequently,  we  obtain 

a  <K2  njVl(q)-l(q)  +  dist{ qN,q)  (  A(q)  _  _  +  Af)  . 

VxV  VxV  \  VxV  ) 

Thus,  given  any  qjV  — >  q  e  Q,  it  follows  that  AN  — >  0  as  N  — >  oo. 

Since  it  can  be  shown  that  the  approximate  solution  A  depends  continuously  on  q, 
solutions  exist  to  the  problem  (ASIDP)^  for  each  N.  Having  established  the  results  of 
Lemma  6,  we  can  now  use  standard  arguments  [10,  11]  to  prove  the  following  theorem. 

Theorem  2:  Suppose  that  the  hypotheses  [H-0]  to  [H-4]  hold  and  let  q^  be  a  solution 
of  the  problem  (ASIDP)W.  Then  the  sequence  jq^  j  admits  a  convergent  subsequence 
qNk  with  qNk  — »  q*  as  TV*,  —■ ►  oo.  Moreover,  q*  is  a  solution  of  the  problem  (SIDP). 

We  now  turn  to  a  particular  implementation  of  the  method  presented  above.  Let  {B^(x i)}^ 
be  the  series  of  B-splinc  functions  ([12]).  To  characterize  the  crack  depth,  the  shape  func¬ 
tion  r(x i)  is  represented  as 

M+ 1 

x2  =  r(x i;q)  =  ^2  qfIB^r(x1)  for  |xi|  <  c  <  di  (59) 

i= o 

with  =  q%+1  =  h0.  In  order  to  ensure  the  compactness  of  the  parameter  set  Q  ([H-0]), 
we  impose  constraints, 

Q  =  {  q  G  Rm  |  /3  <  qi  <h0  i  —  1,  2,  •  •  • ,  M  j  .  (60) 

Remark:  The  defect  shape  function  satisfies  the  hypotheses  [H-l]  to  [H-4], 

For  the  state  approximation,  let  us  choose  U^=]  {'0(V}(li  as  a  set  of  basis  functions  in 
V.  (In  the  calculations  reported  on  below,  we  used  piecewise  linear  finite  elements.) 
That  is,  for  each  N,  are  linearly  independent  and  U^span{ijjf  }]^Ll  is  dense  in 

L2(Gq  U  G\  U  G2s).  Then  the  approximate  subspaces  can  be  chosen  as  HN  =  HN  x  HN 
where  HN  =  span{i/)^}2A] .  Thus  we  can  reconstruct  the  basis  function  by 

f  (ipR ,  0)  for  i  —  1,2,  ■  ■  ■ ,  N 

=\ 

[  (OjV’i-jv)  for  i  =  N  +  1,  N  +  2,  •  •  • ,  2 IV 
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An  approximate  solution  can  be  then  defined  by 


N  2iV 

^4  :=5>W 

i= 1 

where  the  coefficient  vector  wh  =  {uA},2^  is  chosen  such  that,  for  j  —  1,  2,  •  •  •  2N, 

2  N 

«(q  =  L{  q)(^f  )• 

i=l 

Hence  the  system  can  be  approximated  by  solving  the  linear  system 

KN(q)wN  =  fN(  q)  (61) 

where 

[-fC(q)L-  :=  “(q)^,  V’f)  for  b  i  =  i,  2,  •  •  • ,  27V 


and  where 


[fN(q)}j  :=  L(q)(^)  for 

3  =  1,2,-- 

■,27V. 

The  corresponding  output  can  be  computed  as 

zN(q;uJi)  =  HNwN(q;Ui)  1 

(<  =  1,2,-- 

■,M0) 

(62) 

where 

:=  [H\4j  for  i  —  1,2, 

1,2,  •••,21V. 

Our  computational  algorithm  is  to  seek  the  parameter  G  Q  which  minimizes 


1  M0 

FN{ q)  =  -^|wV(q;Wi)  -  zd(ui)  |  .  (63) 

6  Computational  Experiments 

In  the  numerical  experiments,  the  window,  the  suspect  region,  and  the  size  of  conducting 
material  sample  in  Fig.  9  were  preassigned  as 

Window  (domain)  :  d  =  0.350[m]  h  =  0.240[m] 

Suspect  region  :  c  =  0.025 [m] 

Conducting  Material  :  dc  =  0.025 [m]  h0  =  0.020 [m] 

Assuming  that  the  sample  material  was  a  carbon  fiber  reinforced  plastic  (CFRP),  the 
conductivity  ac  and  the  permeability  /i0  were  taken  as  ac  =  1  x  106  [S/m]  and  fi0  = 
47t  x  10~'  [H/m],  respectively.  In  the  numerical  experiments,  we  assume  that  a  single 
frequency  was  used,  i.  e.,,  M0  =  1  in  (57).  The  alternating  current  Js  and  the  applied 
frequency  /  =  u/2tt  were  respectively  given  by  Js  =  40  [mA]  and  /  =  100 [Hz].  The 
number  of  observation  points  was  set  as  m  =  30  and  the  set  of  scanning  positions  of 
SQUID  coil  probe  given  by  (14)  were  chosen  as 

x]  =  —r  +  Oi  ~  e  —  —r  +  Oi  +  e  (i  —  1,  2,  •  •  • ,  30) 
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Figure  10:  The  decomposition  of  sample  material 


where 

Oi  =  0.0125  +  0.025i,  i  =  1,2,  ■••,30, 

and  where  e  denotes  a  sufficiently  small  positive  value  to  restrict  the  effective  area  of 
the  SQUID  pick-up  coil.  The  scanning  line  in  Fig.  9  were  set  as  r  =  0.0375[m]  and 
hf  —  0.003 [m],  respectively. 

To  discretize  the  system  model  by  a  finite-element  method,  the  reference  domain 
was  divided  into  a  finite  number  of  elements  {ek)k=\  and  a  number  Ne(>  Ke) 
of  nodes  defined  by  {xj  =  (x\,  were  selected  in  the  reference  domain.  Each  element 

is  preassigned  as  an  axiprallel  rectangle  with  nodes  at  the  vertices  (see  e.  g.,  [13]).  The 
number  of  finite  elements  and  nodes  in  the  numerical  experiments  reported  in  the  sequel 
were  set  as  Ke  =  3000(=  50x60)  and  Ne  =  3061(=  51x61)  respectively.  Figure  10  depicts 
the  decomposition  of  the  conducting  material  treated  here.  The  number  of  elements  and 
nodes  in  the  material  were  taken  as  200(=  10  x  20)  and  231  (=  11  x  21),  respectively. 
Integration  in  the  element  matrices  KN( q)  and  the  element  vector  fN( q)  were  computed 
numerically  by  a  Gauss-Legendre  formula. 

For  these  test  example  computations,  simulated  data  {zj}  were  generated  by  first 
solving  the  finite- Galer  kin  model  (61)  and  (62).  A  series  of  random  Gaussian  noise  were 
added  to  the  numerical  solutions,  thereby  producing  simulated  noisy  data  for  the  algo¬ 
rithm.  The  essential  difficulties  in  (ASIDP)*  come  from  the  fact  that  real  data  are 
heavily  corrupted  by  observation  noise  and  it  often  occurs  that  the  corresponding  model 
output  data  are  far  from  those  practical  data.  Tikhonov  regularization  is  one  possible 
technique  for  avoiding  these  serious  difficulties  in  computational  efforts.  To  this  end,  we 
adopt  Tikhonov’s  regularization  techniques  to  our  inverse  algorithm.  More  precisely,  a 
regularization  term  is  added  to  the  cost  functional  (63)  (see  [14]  for  more  details): 


q)  =  ^(q)  + 


(64) 
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where 

‘1-1  0  .  0  ' 

0  1  -1  0  •••  0 

[Lc]jj  =  :  0  •••  •••  :  G  i?(M~1)xM 

:  :  0 

.0  0  •••  0  1  -1  . 

In  order  to  select  the  “optimal”  parameter  r],  we  used  L-curve  analysis  (see  [15,  16]).  The 
outline  of  this  approach  in  these  experiments  is  summarized  as  follows: 

[L-curve  analysis]  :  Given  r)  G  [0,  oo),  solve 

F? (qfa) )  =  min  Fn  (<l) >  (65) 

for  each  hxed  rj  G  [0,  oo),  collect  the  points 

(  zN(qv;u)  -  zd( u)  2,  Lcq„  2)  (66) 

and  plot  the  above  points.  Then  the  ”  optimal”  parameter  r/°  can  be  chosen  as  a  point 
that  lies  in  the  ”  corner”  of  the  resulting  curve. 

Evaluation  of  the  vector  gradient  of  the  cost  function  (65)  is  the  computationally  expensive 
part  of  our  algorithm.  This  can  be  accomplished  by  using  a  co-state  approach.  For  the 
numerical  results  reported  in  this  paper,  we  used  the  gradient  projection  method  ([17]) 
which  is  a  particularly  useful  technique  for  optimization  problems  with  linear  inequality 
constraints  such  as  those  in  (60).  The  iterative  algorithm  for  finding  qA'  in  this  experiment 
is  the  same  as  that  in  [5]. 

In  our  test  examples,  we  considered  corrosion  shape  identification  for  conducting  mate¬ 
rial  samples  with  corrosion- like  damage  (of  varying  corrosion  shape)  as  depicted  in  Fig.  6. 
The  parameter  function  (59)  to  be  identified  is  a  piecewise  linear  spline  function.  We 
denote  the  knot  sequence  for  r  by 

-0.025  =  r0M  <  rf  <  •  •  •  <  Tjjf+1  =  0.025 

and  the  unknown  parameter  vector  qAi  =  is  then  given  by 

C  =  r(rtM)  fori  =  1,2 ,  •  •  • ,  M. 

For  the  computational  experiments,  the  dimensions  of  unknown  parameter  vector  q  was 
taken  as  M  =  19.  The  lower  bound  in  (60)  was  chosen  as  (3  —  0.001.  The  initial  guesses 
for  the  parameters  were  given  by 

gf  =  ho(=0.02)  (i  =  1,  2,  •  •  • ,  M) 

which  imply  that  there  exists  no  corrosion  in  the  material.  Carrying  out  a  number  of 
numerical  experiments,  we  checked  the  robustness  of  the  algorithm  with  respect  to  noise 
for  various  kinds  of  corrosion  type  damages.  Here  we  summarize  the  numerical  results  for 
two  typical  examples. 
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True  function 
Initial  Guess 


Figure  11:  True  function  and  estimated  function  in  Example  1  (noise  free) 


Example  1:  In  this  test  example,  we  deal  with  the  shape  recovery  for  the  CFRP  sam¬ 
ple  with  a  symmetric  type  of  corrosion  as  shown  in  Figure  11.  The  values  of  the  true 
parameters  were  given  by 

q\ 9  =  q\g  =  0.0195  g219  =  g}9  =  0.0190  q™  =  q 19  =  0.0185  qf  =  g}9  =  0.0180 
g519  =  q\l  =  0.0175  g619  =  q\\  =  0.0170  q)9  =  q\l  =  0.0165  g819  =  q\%  =  0.0160 
g919  =  q\l  =  0.0155  g|9  =  0.0150. 

The  estimated  parameters  are  listed  in  Table  1  using  data  that  is  with  noise  free,  and  with 
data  containing  5%,  10%, and  20%  relative  noise.  Fig.  12  depicts  the  L-curve  for  5%  noise 
case  and  rf  =  5.0  x  10~12  was  chosen  as  the  optimal  regularization  parameter  in  (64). 
The  successful  results  is  demonstrated  in  Fig.  13  that  depicts  the  estimated  shape  using 
Tikhonov  regularization  with  the  optimal  selection  rf  and  without  the  regularization. 


Example  2:  In  this  example,  we  deal  with  a  somewhat  more  difficult  case  as  compared 
to  that  of  Example-1.  The  type  of  damage  was  given  by  a  non-symmetric  function  as 
illustrated  in  Figure  14.  We  set  the  same  dimension  for  the  unknown  parameter  vector  as 
in  Example  1  and  we  also  used  the  same  knot  sequence.  The  values  of  the  true  parameters 
were  preassigned  in 


ffi19 

=  0.0195 

<?219 

=  0.0190 

%19 

=  0.0185 

=  0.0180 

?519 

%19 

=  0.0170 

$ 

=  0.0165 

%19 

=  0.0160 

ql9 

=  0.0155 

qlo 

0.0175 

0.0150 
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Table  1:  True  value  and  estimated  values  in  Example  1 


9i 

92 

93 

94 

95 

96 

97 

98 

CD 

o 

True  Value  (xlO  2) 

1.95 

1.90 

1.85 

1.80 

1.75 

1.70 

1.65 

1.60 

1.55  1.50  1.55 

Initial  Guess  (xlO-2) 

2.00 

2.00 

2.00 

2.00 

2.00 

2.00 

2.00 

2.00 

2.00  2.00  2.00 

Noise 

iteration  10 

1.97 

1.93 

1.89 

1.84 

1.79 

1.74 

1.69 

1.64 

1.58  1.53  1.58 

Free 

iteration  20 

1.96 

1.92 

1.87 

1.82 

1.77 

1.72 

1.66 

1.61 

1.56  1.48  1.56 

(xlO-2) 

iteration  30 

1.95 

1.91 

1.86 

1.81 

1.76 

1.70 

1.66 

1.61 

1.56  1.48  1.56 

5% 

iteration  10 

2.00 

1.99 

1.96 

1.90 

1.83 

1.76 

1.69 

1.61 

1.57  1.59  1.62 

Noise 

(xlO”2) 

iteration  14 

2.00 

1.99 

1.96 

1.90 

1.83 

1.75 

1.67 

1.59 

1.50  1.48  1.58 

10% 

iteration  8 

2.00 

1.99 

1.96 

1.91 

1.85 

1.78 

1.71 

1.63 

1.63  1.68  1.69 

Noise 

(xl0“2) 

iteration  13 

2.00 

1.99 

1.96 

1.90 

1.83 

1.75 

1.66 

1.57 

1.42  1.53  1.72 

20% 

iteration  10 

2.00 

2.00 

1.99 

1.95 

1.90 

1.84 

1.78 

1.72 

1.74  1.80  1.82 

Noise 

(xl0~2) 

iteration  13 

2.00 

2.00 

1.99 

1.95 

1.89 

1.82 

1.75 

1.67 

1.57  1.70  1.86 

q\2 

9l3 

914 

9l5 

916 

9l7 

918 

919 

True  Value(xl0  2) 

1.60 

1.65 

1.70 

1.75 

1.80 

1.85 

1.90 

1.95 

Initial  Guess(xl0~2) 

2.00 

2.00 

2.00 

2.00 

2.00 

2.00 

2.00 

2.00 

Noise 

iteration  10 

1.64 

1.69 

1.74 

1.79 

1.84 

1.89 

1.93 

1.97 

1.26  x  10"5 

Free 

iteration  20 

1.61 

1.66 

1.72 

1.77 

1.82 

1.87 

1.92 

1.96 

2.62  x  10“6 

(xlO-2) 

iteration  30 

1.61 

1.66 

1.71 

1.76 

1.81 

1.86 

1.91 

1.95 

9.19  x  10~7 

5% 

iteration  10 

1.64 

1.69 

1.77 

1.84 

1.89 

1.94 

1.97 

2.00 

5.11  x  10"5 

Noise 

(xlO’2) 

iteration  14 

1.63 

1.65 

1.73 

1.82 

1.90 

1.94 

1.97 

2.00 

4.16  x  10~5 

10% 

iteration  8 

1.69 

1.73 

1.80 

1.87 

1.92 

1.95 

1.97 

1.99 

9.45  x  10"5 

Noise 

(xl0~2) 

iteration  13 

1.64 

1.62 

1.73 

1.85 

1.92 

1.96 

1.97 

1.99 

7.19  x  10"5 

20% 

iteration  10 

1.81 

1.82 

1.86 

1.91 

1.95 

1.98 

1.99 

1.99 

2.62  x  10~4 

Noise 

(xlO-2) 

iteration  13 

1.81 

1.74 

1.80 

1.88 

1.95 

1.99 

2.00 

1.99 

1.96  x  10-4 
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2.5e-05 


Residual  norm 


Figure  12:  Tikhonov  ”L”  -curve  in  Example  1  (5%  noise) 


True  Function  - 

Estimated  Function  using  Tikhonov  regularization  -  —  - 


Estimated  Function  without  regularization 


Figure  13:  True  function  and  estimated  function  in  Example  1  (5%  noise) 
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True  function 
Initial  Guess 


Figure  14:  True  function  and  estimated  function  in  Example  2  (noise  free) 


qlf  =  0.0145  q\l  =  0.0140  q}g  =  0.0135  q{g  =  0.0130  q}g  =  0.0125 
q\l  =  0.0120  q$  =  0.0140  q}g  =  0.0160  q\9  =  0.0180. 

The  dimension  of  the  unknown  parameter  vector  and  the  number  of  observation  points 
were  chosen  the  same  as  in  Example  1.  Table  2  shows  the  estimated  parameters  for 
the  data,  that  is,  with  noise  free,  and  data  with  5%,  10%,  and  20%  relative  noise.  Fig¬ 
ure  15  depicts  the  Tikhonov  L-curve  for  5%  noise  case.  Estimated  corrosion  shapes  using 
Tikhonov  regularization  and  without  the  regularization  were  demonstrated  in  Fig.  16. 


27 


x2  [m]  Semi-norm 


Figure  15:  Tikhonov  ”L”-curve  in  Example  2  (5%  noise) 


True  Function  - 

Estimated  Function  using  Tikhonov  regularization  —  —  - 


Estimated  Function  without  regularization 


Figure  16:  True  function  and  estimated  function  in  Example  2  (5%  noise) 


Table  2:  True  value  and  estimated  values  in  Example  2 


9i 

92 

93 

94 

95 

96 

97 

98 

CD 

O 

True  Value(xlO  2) 

1.95 

1.90 

1.85 

1.80 

1.75 

1.70 

1.65 

1.60 

1.55  1.50  1.45 

Initial  Guess(x  10”2) 

2.00 

2.00 

2.00 

2.00 

2.00 

2.00 

2.00 

2.00 

2.00  2.00  2.00 

Noise 

iteration  10 

1.99 

1.97 

1.94 

1.90 

1.86 

1.81 

1.77 

1.72 

1.67  1.62  1.57 

Free 

iteration  20 

1.97 

1.94 

1.89 

1.84 

1.79 

1.73 

1.68 

1.63 

1.59  1.55  1.52 

(xlO~2) 

iteration  190 

1.92 

1.93 

1.85 

1.80 

1.74 

1.72 

1.63 

1.58 

1.54  1.49  1.45 

5% 

iteration  9 

2.00 

1.99 

1.99 

1.91 

1.92 

1.80 

1.80 

1.70 

1.66  1.63  1.61 

Noise 

(xlCT2) 

iteration  12 

2.00 

2.00 

1.97 

1.93 

1.87 

1.81 

1.75 

1.68 

1.62  1.57  1.55 

10% 

iteration  10 

2.00 

2.00 

1.96 

1.94 

1.87 

1.83 

1.76 

1.72 

1.61  1.61  1.64 

Noise 

(xl0“2) 

iteration  13 

2.00 

2.00 

1.97 

1.93 

1.87 

1.81 

1.75 

1.68 

1.60  1.50  1.59 

20% 

iteration  10 

2.00 

2.00 

1.97 

1.94 

1.87 

1.83 

1.74 

1.69 

1.52  1.65  1.73 

Noise 

(xl0~2) 

iteration  13 

2.00 

2.00 

1.98 

1.93 

1.87 

1.80 

1.73 

1.63 

1.48  1.49  1.78 

912 

913 

914 

9l5 

916 

9l7 

918 

919 

True  Value(xl0  2) 

1.40 

1.35 

1.30 

1.25 

1.20 

1.40 

1.60 

1.80 

Initial  Guess(xl0~2) 

2.00 

2.00 

2.00 

2.00 

2.00 

2.00 

2.00 

2.00 

Noise 

iteration  10 

1.53 

1.47 

1.26 

1.40 

1.52 

1.60 

1.72 

1.85 

1.76  x  10"4 

Free 

iteration  20 

1.48 

1.38 

1.22 

1.33 

1.52 

1.43 

1.66 

1.80 

7.86  x  10~5 

(xl0~2) 

iteraton  190 

1.42 

1.36 

1.29 

1.19 

1.38 

1.35 

1.55 

1.75 

2.43  x  10"5 

5% 

iteration  9 

1.59 

1.52 

1.37 

1.47 

1.62 

1.70 

1.78 

1.92 

3.13  x  10-4 

Noise 

(xl0~2) 

iteration  12 

1.58 

1.52 

1.31 

1.29 

1.52 

1.72 

1.79 

1.93 

2.24  x  10~4 

10% 

iteration  10 

1.63 

1.51 

1.41 

1.55 

1.70 

1.77 

1.80 

1.92 

4.08  x  10“4 

Noise 

(xlO”2) 

iteration  13 

1.66 

1.54 

1.31 

1.40 

1.65 

1.80 

1.81 

1.91 

3.45  x  10“4 

20% 

iteration  10 

1.66 

1.53 

1.55 

1.68 

1.79 

1.85 

1.85 

1.91 

6.17  x  10~4 

Noise 

(xlO-2) 

iteration  13 

1.73 

1.49 

1.40 

1.58 

1.78 

1.88 

1.85 

1.90 

5.71  x  10-4 
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7  Concluding  Remarks 

A  computational  method  was  considered  for  the  nonlinear  output  least  square  problem 
arising  in  electromagnetic  nondestructive  testing  based  on  SQUID  measurements.  An 
important  feature  of  the  problem  treated  here  was  the  use  of  a  method  of  mapping  tech¬ 
nique  for  identifying  the  corrosion  shape  on  the  back  surface  of  the  conducting  material. 
This  method  was  effectively  used  both  in  the  theoretical  convergence  for  approximating 
identification  problems  and  in  the  related  computational  methods.  The  practical  utility  of 
our  computational  algorithm  is  demonstrated  through  a  series  of  numerical  experiments 
as  summarized  in  the  previous  section.  Carrying  out  a  large  number  of  other  numerical 
tests  in  addition  to  those  reported  for  Examples  1  and  2,  we  report  that  the  algorithm 
performed  quite  well  for  various  noise  levels.  In  the  numerical  experiments,  robustness  of 
the  algorithm  with  respect  to  observation  noise  in  the  simulated  data  was  demonstrated. 
A  graphical  method  was  also  tested  for  determining  the  Tikhonov  parameter  by  using  the 
L-curve  method.  In  all  cases,  our  computational  packages  performed  very  well  up  to  5% 
noise  level.  We  are  currently  pursuing  computational  experiments  with  laboratory  data. 
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